MCMC convergence diagnostics`
Adam Howes
Abstract
Background
Task
out <- readRDS("depends/out.rds")
mcmc <- out$mcmc
bayesplot::color_scheme_set("viridisA")
ggplot2::theme_set(theme_minimal())
#' Rhat
#' Looking for values less than 1.1 here
rhats <- bayesplot::rhat(mcmc)
bayesplot::mcmc_rhat(rhats)
(big_rhats <- rhats[rhats > 1.1])## beta_rho[1] beta_rho[2] beta_alpha[1] beta_alpha[2] logit_phi_rho_x log_sigma_rho_x logit_phi_rho_xs log_sigma_rho_xs logit_phi_rho_a log_sigma_rho_a logit_phi_rho_as log_sigma_rho_as u_rho_x[4]
## 1.436462 2.050365 1.273354 1.462126 1.248576 1.150130 1.494549 1.362138 1.126501 1.151890 1.256854 1.325939 1.171248
## us_rho_x[4] us_rho_x[6] u_rho_xs[1] u_rho_xs[2] u_rho_xs[4] u_rho_xs[11] u_rho_xs[12] u_rho_xs[22] u_rho_xs[27] u_rho_xs[30] us_rho_xs[1] us_rho_xs[2] us_rho_xs[3]
## 1.139724 1.159936 1.108818 1.104644 1.118410 1.113635 1.110012 1.174033 1.201293 1.139680 1.127415 1.138374 1.174318
## us_rho_xs[4] us_rho_xs[5] us_rho_xs[6] us_rho_xs[8] us_rho_xs[9] us_rho_xs[10] us_rho_xs[11] us_rho_xs[12] us_rho_xs[13] us_rho_xs[14] us_rho_xs[15] us_rho_xs[22] us_rho_xs[23]
## 1.133138 1.156225 1.183878 1.133482 1.202853 1.187893 1.199672 1.198588 1.185517 1.125922 1.194364 1.185790 1.110124
## us_rho_xs[24] us_rho_xs[25] us_rho_xs[26] us_rho_xs[27] us_rho_xs[28] us_rho_xs[29] us_rho_xs[30] us_rho_xs[31] us_rho_xs[32] u_rho_a[1] u_rho_a[2] u_rho_a[3] u_rho_a[4]
## 1.238823 1.134849 1.259984 1.296719 1.176014 1.247366 1.324837 1.255172 1.188287 1.429221 1.429893 1.435065 1.430681
## u_rho_a[5] u_rho_a[6] u_rho_a[7] u_rho_a[8] u_rho_a[9] u_rho_a[10] u_rho_as[1] u_rho_as[2] u_rho_as[3] u_rho_as[4] u_rho_as[5] u_rho_as[6] u_rho_as[7]
## 1.432134 1.432472 1.432241 1.436418 1.434241 1.434673 2.017298 2.033510 2.041822 2.038281 2.054323 2.044709 2.034816
## u_rho_as[8] u_rho_as[9] u_rho_as[10] logit_phi_alpha_x log_sigma_alpha_x logit_phi_alpha_xs log_sigma_alpha_xs u_alpha_x[1] u_alpha_x[2] u_alpha_x[3] u_alpha_x[4] u_alpha_x[8] u_alpha_x[9]
## 2.029767 2.036252 1.996995 1.101527 1.410485 1.146635 1.268746 1.263213 1.166815 1.240560 1.177752 1.193106 1.171466
## u_alpha_x[10] u_alpha_x[13] u_alpha_x[14] u_alpha_x[21] u_alpha_x[22] u_alpha_x[23] u_alpha_x[24] u_alpha_x[25] u_alpha_x[26] u_alpha_x[27] u_alpha_x[28] u_alpha_x[30] us_alpha_x[1]
## 1.187447 1.273750 1.117343 1.432030 1.195067 1.190772 1.272183 1.363360 1.142784 1.178192 1.109855 1.130643 1.168862
## us_alpha_x[2] us_alpha_x[3] us_alpha_x[4] us_alpha_x[5] us_alpha_x[6] us_alpha_x[8] us_alpha_x[12] us_alpha_x[15] us_alpha_x[21] us_alpha_x[22] us_alpha_x[24] us_alpha_x[25] us_alpha_x[29]
## 1.133704 1.197043 1.264392 1.112401 1.198646 1.100455 1.118491 1.118520 1.189470 1.245977 1.164936 1.132086 1.162871
## us_alpha_x[30] u_alpha_xs[1] u_alpha_xs[8] u_alpha_xs[11] u_alpha_xs[13] u_alpha_xs[16] u_alpha_xs[17] u_alpha_xs[18] u_alpha_xs[19] u_alpha_xs[21] u_alpha_xs[22] u_alpha_xs[23] u_alpha_xs[24]
## 1.136288 1.115010 1.139084 1.103883 1.162435 1.145249 1.133684 1.299266 1.209079 1.405892 1.201692 1.144698 1.275732
## u_alpha_xs[25] u_alpha_xs[28] u_alpha_xs[29] u_alpha_xs[32] us_alpha_xs[1] us_alpha_xs[2] us_alpha_xs[3] us_alpha_xs[4] us_alpha_xs[5] us_alpha_xs[6] us_alpha_xs[8] us_alpha_xs[9] us_alpha_xs[10]
## 1.257591 1.187026 1.390206 1.111991 1.158874 1.175500 1.192279 1.138272 1.181611 1.150340 1.217743 1.187316 1.192677
## us_alpha_xs[11] us_alpha_xs[12] us_alpha_xs[13] us_alpha_xs[15] us_alpha_xs[17] us_alpha_xs[18] us_alpha_xs[19] us_alpha_xs[20] us_alpha_xs[21] us_alpha_xs[22] us_alpha_xs[23] us_alpha_xs[24] us_alpha_xs[25]
## 1.180764 1.122971 1.176139 1.124448 1.145104 1.138565 1.131682 1.104125 1.154424 1.203935 1.146553 1.133402 1.145399
## us_alpha_xs[26] u_alpha_a[1] u_alpha_a[2] u_alpha_a[3] u_alpha_a[4] u_alpha_a[5] u_alpha_a[6] u_alpha_a[7] u_alpha_a[8] u_alpha_a[9] u_alpha_a[10] u_alpha_a[11] u_alpha_a[12]
## 1.115193 1.237430 1.260247 1.265380 1.224632 1.234161 1.248809 1.257960 1.259736 1.262103 1.253853 1.255800 1.234374
## u_alpha_a[13] u_alpha_as[1] u_alpha_as[2] u_alpha_as[3] u_alpha_as[4] u_alpha_as[5] u_alpha_as[6] u_alpha_as[7] u_alpha_as[8] u_alpha_as[9] u_alpha_as[10] u_alpha_xa[13] u_alpha_xa[21]
## 1.218461 1.405071 1.431038 1.426739 1.435003 1.462654 1.469958 1.488702 1.483664 1.495276 1.430269 1.107958 1.335780
## u_alpha_xa[24] u_alpha_xa[25] u_alpha_xa[29] log_sigma_lambda_x ui_lambda_x[1] ui_lambda_x[10] ui_lambda_x[11] ui_lambda_x[13] ui_lambda_x[14] ui_lambda_x[22] ui_lambda_x[23] ui_lambda_x[32] log_sigma_ancalpha_x
## 1.259357 1.108961 1.204265 2.309844 1.212403 1.427155 1.258550 1.187641 1.158190 1.139450 1.131882 1.146537 1.470755
## ui_anc_rho_x[15] ui_anc_rho_x[25] ui_anc_rho_x[29] ui_anc_alpha_x[16] ui_anc_alpha_x[18] ui_anc_alpha_x[19] ui_anc_alpha_x[24] ui_anc_alpha_x[25] ui_anc_alpha_x[29] ui_anc_alpha_x[32] log_or_gamma[3] log_or_gamma[5] log_or_gamma[8]
## 1.101091 1.172969 1.114979 1.195326 1.303179 1.144301 1.125907 1.151146 1.269189 1.129052 1.124149 1.113333 1.101212
## log_or_gamma[9] log_or_gamma[10] log_or_gamma[12] log_or_gamma[18] log_or_gamma[20] log_or_gamma[22] log_or_gamma[31] lp__
## 1.123776 1.127447 1.141006 1.102396 1.114067 1.128760 1.113397 1.567209length(big_rhats) / length(rhats)## [1] 0.4126016#' ESS ratio
#' Worry about values less than 0.1 here... which they all are :clown:
ratios <- bayesplot::neff_ratio(mcmc)
bayesplot::mcmc_neff(ratios)
#' Autocorrelation
#' High values of autocorrelation in the chains
bayesplot::mcmc_acf(mcmc, pars = vars(starts_with("beta")))
#' Univariate traceplots
#' Assess these visually
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("beta")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("logit")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("log_sigma")))
#' Prevalence model
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_xs[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_rho_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_rho_xs[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_a[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_as[")))
#' ART model
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_xs[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_alpha_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_alpha_xs[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_a[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_as[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_xa[")))
#' Other
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_lambda_x[")))
#' ANC model
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_anc_rho_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_anc_alpha_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("log_or_gamma["))) #' N.B. these are from the ANC attendance model
#' Pairs plots
#' Prior suspicion (from Jeff, Tim, Rachel) that the ART attendance model is unidentifiable
#' Let's have a look at the pairs plot for neighbouring districts and the log_or_gamma parameter
nb <- area_merged %>%
filter(area_level == max(area_level)) %>%
bsae::sf_to_nb()
neighbours_pairs_plot <- function(par, i) {
neighbour_pars <- paste0(par, "[", c(i, nb[[i]]), "]")
bayesplot::mcmc_pairs(mcmc, pars = neighbour_pars, diag_fun = "hist", off_diag_fun = "hex")
}
area_merged %>%
filter(area_level == max(area_level)) %>%
print(n = Inf)## Simple feature collection with 32 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 32.67152 ymin: -17.12721 xmax: 35.91505 ymax: -9.364003
## Geodetic CRS: WGS 84
## # A tibble: 32 × 12
## area_id area_name area_level parent_area_id spectrum_region_code area_sort_order center_x center_y area_level_label display naomi_level geometry
## * <chr> <chr> <int> <chr> <dbl> <int> <dbl> <dbl> <chr> <lgl> <lgl> <MULTIPOLYGON [°]>
## 1 MWI_4_1_demo Chitipa 4 MWI_3_1_demo 0 38 33.4 -9.98 District + Metro TRUE TRUE (((33.63988 -9.603465, 33.64077 -9.611695, 33.64722 -9.620278, 33.65034 -9.630591, 3...
## 2 MWI_4_2_demo Karonga 4 MWI_3_2_demo 0 39 33.9 -10.1 District + Metro TRUE TRUE (((34.17391 -10.58124, 34.16676 -10.57488, 34.15757 -10.57275, 34.15369 -10.56901, 3...
## 3 MWI_4_3_demo Rumphi 4 MWI_3_3_demo 0 40 33.8 -10.8 District + Metro TRUE TRUE (((33.38671 -11.14385, 33.37861 -11.14527, 33.37376 -11.13834, 33.37219 -11.12912, 3...
## 4 MWI_4_4_demo Mzuzu City 4 MWI_3_4_demo 0 41 34.0 -11.4 District + Metro TRUE TRUE (((33.96686 -11.49152, 33.96384 -11.48679, 33.95983 -11.483, 33.95649 -11.47607, 33....
## 5 MWI_4_5_demo Nkhata Bay 4 MWI_3_5_demo 0 42 34.0 -11.7 District + Metro TRUE TRUE (((34.01984 -12.23471, 34.01532 -12.23359, 34.00829 -12.2278, 34.0042 -12.22929, 33....
## 6 MWI_4_6_demo Mzimba 4 MWI_3_4_demo 0 43 33.6 -11.8 District + Metro TRUE TRUE (((33.53559 -12.37776, 33.54323 -12.37243, 33.54602 -12.36943, 33.54572 -12.36572, 3...
## 7 MWI_4_7_demo Likoma 4 MWI_3_6_demo 0 44 34.7 -12.1 District + Metro TRUE TRUE (((34.73153 -12.03097, 34.73347 -12.03236, 34.73292 -12.03681, 34.73486 -12.03875, 3...
## 8 MWI_4_8_demo Nkhotakota 4 MWI_3_7_demo 0 45 34.1 -12.8 District + Metro TRUE TRUE (((34.12471 -12.3983, 34.13183 -12.39831, 34.14007 -12.3996, 34.14636 -12.40232, 34....
## 9 MWI_4_9_demo Kasungu 4 MWI_3_8_demo 0 46 33.5 -13.0 District + Metro TRUE TRUE (((33.32294 -13.61568, 33.31926 -13.61632, 33.31408 -13.6193, 33.3055 -13.61247, 33....
## 10 MWI_4_10_demo Ntchisi 4 MWI_3_9_demo 0 47 33.9 -13.3 District + Metro TRUE TRUE (((34.1272 -13.50292, 34.12042 -13.50349, 34.11702 -13.50761, 34.1156 -13.51225, 34....
## 11 MWI_4_11_demo Dowa 4 MWI_3_10_demo 0 48 33.7 -13.5 District + Metro TRUE TRUE (((33.83662 -13.7742, 33.83127 -13.77034, 33.82476 -13.75987, 33.82003 -13.75745, 33...
## 12 MWI_4_12_demo Salima 4 MWI_3_11_demo 0 49 34.4 -13.8 District + Metro TRUE TRUE (((34.55063 -14.10458, 34.54295 -14.10673, 34.53391 -14.1058, 34.52168 -14.10958, 34...
## 13 MWI_4_13_demo Mchinji 4 MWI_3_12_demo 0 50 33.0 -13.8 District + Metro TRUE TRUE (((32.92279 -13.37107, 32.92497 -13.36179, 32.92436 -13.3568, 32.94433 -13.35478, 32...
## 14 MWI_4_14_demo Lilongwe City 4 MWI_3_13_demo 0 51 33.8 -13.9 District + Metro TRUE TRUE (((33.75722 -13.75277, 33.76022 -13.75845, 33.77182 -13.76671, 33.77989 -13.76876, 3...
## 15 MWI_4_15_demo Lilongwe 4 MWI_3_13_demo 0 52 33.5 -14.0 District + Metro TRUE TRUE (((33.70239 -14.54722, 33.71319 -14.57181, 33.70703 -14.57638, 33.70465 -14.58096, 3...
## 16 MWI_4_16_demo Dedza 4 MWI_3_14_demo 0 53 34.3 -14.2 District + Metro TRUE TRUE (((34.40013 -14.412, 34.39294 -14.40303, 34.39324 -14.3956, 34.38913 -14.39334, 34.3...
## 17 MWI_4_17_demo Ntcheu 4 MWI_3_15_demo 0 54 34.8 -14.8 District + Metro TRUE TRUE (((34.80996 -15.32603, 34.80709 -15.32502, 34.8033 -15.3176, 34.80719 -15.31364, 34....
## 18 MWI_4_18_demo Mangochi 4 MWI_3_16_demo 0 55 35.3 -14.1 District + Metro TRUE TRUE (((34.97112 -14.77456, 34.96737 -14.7704, 34.95998 -14.75692, 34.95798 -14.75167, 34...
## 19 MWI_4_19_demo Machinga 4 MWI_3_17_demo 0 56 35.6 -14.9 District + Metro TRUE TRUE (((35.16126 -15.20686, 35.16416 -15.19912, 35.16824 -15.19152, 35.1692 -15.18516, 35...
## 20 MWI_4_20_demo Balaka 4 MWI_3_18_demo 0 57 35.1 -15.0 District + Metro TRUE TRUE (((35.03931 -15.3229, 35.00434 -15.32357, 34.98965 -15.3241, 34.91631 -15.32533, 34....
## 21 MWI_4_21_demo Zomba City 4 MWI_3_19_demo 0 58 35.3 -15.4 District + Metro TRUE TRUE (((35.34807 -15.41108, 35.34317 -15.4197, 35.3424 -15.42337, 35.33797 -15.42816, 35....
## 22 MWI_4_22_demo Zomba 4 MWI_3_19_demo 0 59 35.4 -15.4 District + Metro TRUE TRUE (((35.84263 -15.45047, 35.83955 -15.44715, 35.83957 -15.44162, 35.83749 -15.43646, 3...
## 23 MWI_4_23_demo Phalombe 4 MWI_3_20_demo 0 60 35.7 -15.7 District + Metro TRUE TRUE (((35.80985 -15.93286, 35.80041 -15.93066, 35.79598 -15.93278, 35.78529 -15.92961, 3...
## 24 MWI_4_24_demo Mulanje 4 MWI_3_21_demo 0 61 35.5 -15.9 District + Metro TRUE TRUE (((35.30769 -16.21557, 35.30645 -16.21306, 35.29699 -16.20214, 35.28732 -16.19913, 3...
## 25 MWI_4_25_demo Neno 4 MWI_3_22_demo 0 62 34.7 -15.5 District + Metro TRUE TRUE (((34.73322 -15.81642, 34.72392 -15.81112, 34.71292 -15.81164, 34.71111 -15.80616, 3...
## 26 MWI_4_26_demo Blantyre 4 MWI_3_23_demo 0 63 34.9 -15.7 District + Metro TRUE TRUE (((34.90176 -16.01093, 34.8995 -16.01271, 34.89635 -16.00977, 34.8893 -16.01434, 34....
## 27 MWI_4_27_demo Mwanza 4 MWI_3_24_demo 0 64 34.5 -15.6 District + Metro TRUE TRUE (((34.73546 -15.81671, 34.73438 -15.81708, 34.52702 -15.81355, 34.52308 -15.80903, 3...
## 28 MWI_4_28_demo Chiradzulu 4 MWI_3_25_demo 0 65 35.2 -15.8 District + Metro TRUE TRUE (((35.31332 -15.98932, 35.30383 -16.00106, 35.30011 -16.00326, 35.2948 -16.00073, 35...
## 29 MWI_4_29_demo Blantyre City 4 MWI_3_23_demo 0 66 35.1 -15.8 District + Metro TRUE TRUE (((35.10285 -15.85812, 35.1012 -15.86633, 35.09855 -15.87003, 35.09046 -15.87081, 35...
## 30 MWI_4_30_demo Thyolo 4 MWI_3_26_demo 0 67 35.1 -16.1 District + Metro TRUE TRUE (((35.26228 -16.40585, 35.25777 -16.4048, 35.25433 -16.40626, 35.24886 -16.4003, 35....
## 31 MWI_4_31_demo Chikwawa 4 MWI_3_27_demo 0 68 34.7 -16.2 District + Metro TRUE TRUE (((34.9105 -16.69166, 34.90701 -16.68995, 34.88869 -16.68457, 34.8867 -16.6869, 34.8...
## 32 MWI_4_32_demo Nsanje 4 MWI_3_28_demo 0 69 35.1 -16.7 District + Metro TRUE TRUE (((34.9105 -16.69166, 34.91037 -16.688, 34.90659 -16.68455, 34.90675 -16.67798, 34.9...neighbours_pairs_plot("log_or_gamma", 5) #' Nkhata Bay and neighbours
neighbours_pairs_plot("log_or_gamma", 26) #' Blantyre and neighbours
#' NUTS specific assessment
np <- bayesplot::nuts_params(mcmc)
#' Number of divergent transitions
#' In this instance it's zero, so no need to do further divergent transition investigation
np %>%
filter(Parameter == "divergent__") %>%
summarise(n_divergent = sum(Value))bayesplot::mcmc_nuts_divergence(np, bayesplot::log_posterior(mcmc))
#' Energy plots
#' Ideally these two histograms would be the same (Betancourt 2017)
#' The histograms are quite different, in a way suggesting the chains may not be
#' fully exploring the tails of the target distribution
bayesplot::mcmc_nuts_energy(np)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.